Morphology Module

History

current version 1.2 - 22nd April 2024

version date comment
0.1 17/Aug/2009 Original code
0.2 23/Oct/2010 Compute flow accumulation grid
0.3 04/Dec/2012 Compute longest flow length
0.4 03/Apr/2013 Derive aspect subroutine
0.5 10/Apr/2013 Derive curvature subroutine
0.6 13/Apr/2013 Fixed DeriveSlope to match ArcGis algorithm
0.7 16/Apr/2013 DrainageDensity rooutine
0.8 03/May/2013 TopographicIndex routine
0.9 10/Jun/2013 CurvatureMicroMet routine
1.0 09/Oct/2013 LongestFlowPathSlope
1.1 03/Feb/2021 Compute Centroid
1.2 22/Apr/2024 ESRI and GRASS flow direction can be used

License

license: GNU GPL http://www.gnu.org/licenses/

This file is part of

MOSAICO -- MOdular library for raSter bAsed hydrologIcal appliCatiOn.

Copyright (C) 2011 Giovanni Ravazzani

Code Description

Language: Fortran 90.

Software Standards: "European Standards for Writing and
Documenting Exchangeable Fortran 90 Code".

Module Description

library to deal with river and basin morphology


Uses


Variables

Type Visibility Attributes Name Initial
integer(kind=short), public :: E

East

integer(kind=short), public, parameter :: ESRI = 1
integer(kind=short), public, parameter :: GRASS = 2
integer(kind=short), public :: N

North

integer(kind=short), public :: NE

North-East
integer(kind=short), public :: NW

North-West
integer(kind=short), public :: S

South

integer(kind=short), public :: SE

South-East
integer(kind=short), public :: SW

South-West
integer(kind=short), public :: W

West

integer(kind=short), public :: flowDirectionConvention

Interfaces

public interface Centroid

  • private subroutine CentroidGridReal(grid, ccoord)

    returns the coordinate of the centroid given a grid_integer mask, in the same coordinate reference system of grid_integer

    Arguments

    Type IntentOptional Attributes Name
    type(grid_real), intent(in) :: grid

    input grid
    type(Coordinate), intent(out) :: ccoord

    coordinate of computed centroid

  • private subroutine CentroidGridInteger(grid, ccoord)

    returns the coordinate of the centroid given a grid_integer mask, in the same coordinate reference system of grid_integer

    Arguments

    Type IntentOptional Attributes Name
    type(grid_integer), intent(in) :: grid

    input grid
    type(Coordinate), intent(out) :: ccoord

    coordinate of computed centroid

Functions

public function CellIsSpring(row, col, flowDir) result(spring)

find the cells that are springs, defined as those cells that have not any other cells upstream

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: row
integer, intent(in) :: col
type(grid_integer), intent(in) :: flowDir

Return Value logical

public function CheckOutlet(icurrent, jcurrent, iDown, jDown, flowDirection)

if the downstream cell is a nodata value or out of grid space, or points toward the current cell, the current cell is the basin outlet.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: icurrent
integer, intent(in) :: jcurrent
integer, intent(in) :: iDown
integer, intent(in) :: jDown
type(grid_integer), intent(in) :: flowDirection

Return Value logical

public function DrainageDensity(channel, fdir, x, y, flowAcc) result(dd)

Drainage Density (1/m) of a basin is the total line length of all streams divided by basin area. If flow accumulation is not given, it is computed from flow direction. If coordinate of closing section is not given drainage density is computed for the entire grid.

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(in) :: channel

mask of channel cells. NoData for non channel cells
type(grid_integer), intent(in) :: fdir

flow direction
real(kind=float), intent(in), optional :: x

coordinate of basin closing section

real(kind=float), intent(in), optional :: y

coordinate of basin closing section

type(grid_real), intent(in), optional :: flowAcc

flow accumulation

Return Value real(kind=float)

public function LongestFlowLength(fdir, x, y) result(lmax)

compute longest flow length (m)

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(in) :: fdir
real(kind=float), intent(in) :: x
real(kind=float), intent(in) :: y

Return Value real(kind=float)

public function LongestFlowPathSlope(fdir, dem, x, y, outsection, headsection) result(slope)

compute slope of longest flow path (m/m)

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(in) :: fdir
type(grid_real), intent(in) :: dem
real(kind=float), intent(in) :: x
real(kind=float), intent(in) :: y
type(Coordinate), intent(out), optional :: outsection

coordinate of outlet section
type(Coordinate), intent(out), optional :: headsection

coordinate of head section

Return Value real(kind=float)

public function Perimeter(grid) result(p)

find the cells that are springs, defined as those cells that have not any other cells upstream

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(in) :: grid

Return Value real(kind=float)


Subroutines

public subroutine BasinDelineate(fdir, x, y, mask)

compute mask of river basin given map of flow direction and the coordinate of the outlet point

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(in) :: fdir
real(kind=float), intent(in) :: x

coordinate of outlet
real(kind=float), intent(in) :: y

coordinate of outlet
type(grid_integer), intent(inout) :: mask

public subroutine CurvatureMicroMet(dem, curve_len_scale, curvature)

Calculate curvature according to method implemented in MICROMET

Read more…

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(in) :: dem
real(kind=float), intent(in) :: curve_len_scale
type(grid_real), intent(out) :: curvature

public subroutine DeriveAspect(dem, aspect)

Aspect identifies the downslope direction of the maximum rate of change in value from each cell to its neighbors. Aspect can be thought of as the slope direction. The values of the output raster will be the compass direction of the aspect. A moving window visits each cell in the input raster and for each cell in the center of the window, an aspect value is calculated using an algorithm that incorporates the values of the cell's eight neighbors. If any neighborhood cells are NoData, they are first assigned the value of the center cell, then the aspect is computed. The cells are identified as letters 'a' to 'i', with 'e' representing the cell for which the aspect is being calculated. The rate of change in the x direction for cell 'e' is calculated with the following algorithm:

Read more…

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(in) :: dem
type(grid_real), intent(out) :: aspect

public subroutine DeriveCurvature(dem, curvature, profileCurvature, planCurvature)

Calculates the curvature [1/m] of a raster surface, optionally including profile and plan curvature. If any neighborhood cells are NoData, they are first assigned the value of the center cell. It calculates the second derivative value of the input surface on a cell-by-cell basis. For each cell, a fourth-order polynomial of the form:

Read more…

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(in) :: dem
type(grid_real), intent(out) :: curvature
type(grid_real), intent(out), optional :: profileCurvature
type(grid_real), intent(out), optional :: planCurvature

public subroutine DeriveSlope(dem, slope)

For each cell, the routine calculates the maximum rate of change in value from that cell to its neighbors. Basically, the maximum change in elevation over the distance between the cell and its eight neighbors identifies the steepest downhill descent from the cell. Conceptually, the routine fits a plane to the z-values of a 3 x 3 cell neighborhood around the processing or center cell. The slope value of this plane is calculated using the average maximum technique. The direction the plane faces is the aspect for the processing cell. If there is a cell location in the neighborhood with a NoData z-value, the z-value of the center cell will be assigned to the location.

Read more…

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(in) :: dem
type(grid_real), intent(out) :: slope

public subroutine DownstreamCell(iin, jin, dir, is, js, dx, grid)

returns the position (is,js) of the downstream cell and, optionally, the flow path length, considering cardinal and diagonal direction

Arguments

Type IntentOptional Attributes Name
integer(kind=short), intent(in) :: iin

current cell
integer(kind=short), intent(in) :: jin

current cell
integer(kind=long), intent(in) :: dir

flow direction

integer(kind=short), intent(out) :: is

downstream cell
integer(kind=short), intent(out) :: js

downstream cell
real(kind=float), intent(out), optional :: dx

flow path length [m]
type(grid_integer), intent(in), optional :: grid

used to define coordinate reference system

public subroutine FlowAccumulation(fdir, facc)

compute map of flow accumulation (m2) Input grid: flow direction

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(in) :: fdir
type(grid_real), intent(out) :: facc

public subroutine HortonOrders(flowDirection, orders, basinOrder)

returns a grid_integer containing Horton orders. Horton orders are computed on the entire space-filled basin.

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(in) :: flowDirection
type(grid_integer), intent(inout) :: orders
integer, intent(out), optional :: basinOrder

the maximum order of the basin

public subroutine SetFlowDirectionConvention(string)

Set flow direction convention:

Read more…

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: string

public subroutine TopographicIndex(dem, fdir, facc, xcoord, ycoord, tivalue, tigrid)

Topographic index TI takes into account both a local slope geometry and site location in the landscape combining data on slope steepness S and specific catchment area As:

Read more…

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(in) :: dem

digital elevation model
type(grid_integer), intent(in) :: fdir

flow direction
type(grid_real), intent(in), optional :: facc

flow accumulation in cell unit
real(kind=float), intent(in), optional :: xcoord

x coordinate of cell for computing index
real(kind=float), intent(in), optional :: ycoord

y coordinate of cell for computing index
real(kind=float), intent(out), optional :: tivalue

index of cell with coordinates (x,y)
type(grid_real), intent(out), optional :: tigrid

topographic index

public subroutine UpstreamCell(i, j, flowdir, flowacc, is, js, found, dx)

returns the position (is,js) of the upstream cell and, optionally, the flow path length, considering cardinal and diagonal direction

Arguments

Type IntentOptional Attributes Name
integer(kind=short), intent(in) :: i

current cell
integer(kind=short), intent(in) :: j

current cell
type(grid_integer), intent(in) :: flowdir

flow direction map
type(grid_real), intent(in) :: flowacc

flow accumulation map
integer(kind=short), intent(out) :: is

upstream cell
integer(kind=short), intent(out) :: js

upstream cell
logical, intent(out) :: found

upsteam cell found
real(kind=float), intent(out), optional :: dx

flow path length [m]

private recursive subroutine BasinArea(r, c, fdir, area)

compute basin area (m2)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: r
integer, intent(in) :: c
type(grid_integer), intent(in) :: fdir
real(kind=float), intent(inout) :: area

private recursive subroutine BasinMask(basin, fdir, r, c)

search for cells included in the river basin

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(inout) :: basin
type(grid_integer), intent(in) :: fdir
integer, intent(in) :: r
integer, intent(in) :: c

private subroutine CentroidGridInteger(grid, ccoord)

returns the coordinate of the centroid given a grid_integer mask, in the same coordinate reference system of grid_integer

Arguments

Type IntentOptional Attributes Name
type(grid_integer), intent(in) :: grid

input grid
type(Coordinate), intent(out) :: ccoord

coordinate of computed centroid

private subroutine CentroidGridReal(grid, ccoord)

returns the coordinate of the centroid given a grid_integer mask, in the same coordinate reference system of grid_integer

Arguments

Type IntentOptional Attributes Name
type(grid_real), intent(in) :: grid

input grid
type(Coordinate), intent(out) :: ccoord

coordinate of computed centroid

private subroutine ConfluenceIsAround(row, col, i, j, flowDir, confluence, orders, Norder)

Scan the eigth cells surrounding the center cell (row,col) (neglecting the upstream cell (i,j) ) to find if a confluence is present. The confluence cell must be of the same order or not yet defined.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: row
integer, intent(in) :: col
integer, intent(in) :: i
integer, intent(in) :: j
type(grid_integer), intent(in) :: flowDir
logical, intent(out) :: confluence
type(grid_integer), intent(in) :: orders
integer, intent(in) :: Norder